{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "# Running the Direct Fidelity Estimation (DFE) algorithm\n",
    "This example walks through the steps of running the direct fidelity estimation (DFE) algorithm as described in these two papers:\n",
    "\n",
    "* Direct Fidelity Estimation from Few Pauli Measurements (https://arxiv.org/abs/1104.4695)\n",
    "* Practical characterization of quantum devices without tomography (https://arxiv.org/abs/1104.3835)\n",
    "\n",
    "Optimizations for Clifford circuits are based on a tableau-based simulator:\n",
    "* Improved Simulation of Stabilizer Circuits (https://arxiv.org/pdf/quant-ph/0406196.pdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "bd9529db1c0b"
   },
   "outputs": [],
   "source": [
    "try:\n",
    "    import cirq\n",
    "except ImportError:\n",
    "    print(\"installing cirq...\")\n",
    "    !pip install --quiet cirq\n",
    "    print(\"installed cirq.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import Cirq, DFE, and create a circuit\n",
    "import cirq\n",
    "from cirq.contrib.svg import SVGCircuit\n",
    "import examples.direct_fidelity_estimation as dfe\n",
    "\n",
    "qubits = cirq.LineQubit.range(3)\n",
    "circuit = cirq.Circuit(cirq.CNOT(qubits[0], qubits[2]),\n",
    "                       cirq.Z(qubits[0]),\n",
    "                       cirq.H(qubits[2]),\n",
    "                       cirq.CNOT(qubits[2], qubits[1]))\n",
    "\n",
    "SVGCircuit(circuit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We then create a sampler. For this example, we use a simulator but the code can accept a hardware sampler.\n",
    "noise = cirq.ConstantQubitNoiseModel(cirq.depolarize(0.1))\n",
    "sampler = cirq.DensityMatrixSimulator(noise=noise)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We run the DFE:\n",
    "estimated_fidelity, intermediate_results = dfe.direct_fidelity_estimation(\n",
    "        circuit,\n",
    "        qubits,\n",
    "        sampler,\n",
    "        n_measured_operators=None,  # None=returns all the Pauli strings\n",
    "        samples_per_term=0)  # 0=use dense matrix simulator\n",
    "\n",
    "print('Estimated fidelity: %.2f' % (estimated_fidelity))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# What is happening under the hood?\n",
    "Now, let's look at the `intermediate_results` and correlate what is happening in the code with the papers. The definition of fidelity is:\n",
    "$$\n",
    "F = F(\\hat{\\rho},\\hat{\\sigma}) = \\mathrm{Tr} \\left(\\hat{\\rho} \\hat{\\sigma}\\right)\n",
    "$$\n",
    "where $\\hat{\\rho}$ is the theoretical pure state and $\\hat{\\sigma}$ is the actual state. The idea of DFE is to write fidelity as:\n",
    "$$F= \\sum _i \\frac{\\rho _i \\sigma _i}{d}$$\n",
    "\n",
    "where $d=4^{\\mathit{number-of-qubits}}$, $\\rho _i = \\mathrm{Tr} \\left( \\hat{\\rho} P_i \\right)$, and $\\sigma _i = \\mathrm{Tr} \\left(\\hat{\\sigma} P_i \\right)$. Each of the $P_i$ is a Pauli operator. We can then finally rewrite the fidelity as:\n",
    "\n",
    "$$F= \\sum _i Pr(i) \\frac{\\sigma _i}{\\rho_i}$$\n",
    "\n",
    "with $Pr(i) = \\frac{\\rho_i ^2}{d}$, which is a probability-like set of numbers (between 0.0 and 1.0 and they add up to 1.0).\n",
    "\n",
    "One important question is how do we choose these Pauli operators $P_i$? It depends on whether the circuit is Clifford or not. In case it is, we know that there are \"only\" $2^{\\mathit{number-of-qubits}}$ operators for which $Pr(i)$ is non-zero. In fact, we know that they are all equiprobable with $Pr(i) = \\frac{1}{2^{\\mathit{number-of-qubits}}}$. The code does detect the Cliffordness automatically and switches to this mode. In case the circuit is not Clifford, the code just uses all the operators.\n",
    "\n",
    "Let's inspect that in the case of our example, we do see the Pauli operators with equiprobability (i.e. the $\\rho_i$):\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for pauli_trace in intermediate_results.pauli_traces:\n",
    "    print('Probability %.3f\\tPauli: %s' % (pauli_trace.Pr_i, pauli_trace.P_i))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Yay! We do see 8 entries (we have 3 qubits) with all the same 1/8 probability. What if we had a 23 qubit circuit? In this case, that would be quite many of them. That is where the parameter `n_measured_operators` becomes useful. If it is set to `None` we return *all* the Pauli strings (regardless of whether the circuit is Clifford or not). If set to an integer, we randomly sample the Pauli strings.\n",
    "\n",
    "Then, let's actually look at the measurements, i.e. $\\sigma_i$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for trial_result in intermediate_results.trial_results:\n",
    "    print('rho_i=%.3f\\tsigma_i=%.3f\\tPauli:%s' % (trial_result.pauli_trace.rho_i, trial_result.sigma_i, trial_result.pauli_trace.P_i))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How are these measurements chosen? Since we had set `n_measured_operators=None`, all the measurements are used. If we had set the parameter to an integer, we would only have a subset to start from. We would then, as per the algorithm, sample from this set with replacement according to the probability distribution of $Pr(i)$ (for Clifford circuits, the probabilities are all the same, but for non-Clifford circuits, it means we favor more probable Pauli strings)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What about the parameter `samples_per_term`? Remember that the code can handle both a sampler or use a simulator. If we use a sampler, then we can repeat the measurements `samples_per_term` times. In our case, we use a dense matrix simulator and thus we keep that parameter set to `0`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How do we bound the variance of the fidelity when the circuit is Clifford?\n",
    "Recall that the formula for DFE is:\n",
    "$$F= \\sum _i Pr(i) \\frac{\\sigma _i}{\\rho_i}$$\n",
    "\n",
    "But for Clifford circuits, we have $Pr(i) = \\frac{1}{d}$ and $\\rho_i = 1$ and thus the formula becomes:\n",
    "$$F= \\frac{1}{d} \\sum _i \\sigma _i$$\n",
    "\n",
    "If we estimate by randomly sampling $N$ values for the indicies $i$ for $\\sigma_i$ we get:\n",
    "$$\\hat{F} = \\frac{1}{N} \\sum_{j=1}^N \\sigma _{i(j)}$$\n",
    "\n",
    "Using the Bhatia–Davis inequality ([A Better Bound on the Variance, Rajendra Bhatia and Chandler Davis](https://www.jstor.org/stable/2589180)) and the fact that $0 \\le \\sigma_i \\le 1$, we have the variance of:\n",
    "$$\\mathrm{Var}\\left[ \\hat{F} \\right] \\le \\frac{(1 - F)F}{N}$$\n",
    "\n",
    "$$\\mathrm{StdDev}\\left[ \\hat{F} \\right] \\le \\sqrt{\\frac{(1 - F)F}{N}}$$\n",
    "\n",
    "In particular, since $0 \\le F \\le 1$ we have:\n",
    "$$\\mathrm{StdDev}\\left[ \\hat{F} \\right] \\le \\sqrt{\\frac{(1 - \\frac{1}{2})\\frac{1}{2}}{N}}$$\n",
    "\n",
    "$$\\mathrm{StdDev}\\left[ \\hat{F} \\right] \\le \\frac{1}{2 \\sqrt{N}}$$"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.4rc1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
